{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import string\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = \"/Users/aavattikutis/Documents/epidemicmodel/cccruns/fits/fit2/tables/\"\n",
    "models = [\"fulllinearmodel_fit_table.csv\",\"reducedlinearmodelNegBinom_fit_table.csv\",\n",
    "          \"reducedlinearmodelq0_fit_table.csv\",\"reducedlinearmodelq0ctime_fit_table.csv\",\n",
    "         \"nonlinearmodelq0ctime_fit_table.csv\"]\n",
    "\n",
    "\n",
    "# model_hash = {}\n",
    "# k = -1\n",
    "# for model in models:\n",
    "#     k += 1\n",
    "#     model_hash[model] = string.ascii_uppercase[k]\n",
    "\n",
    "# df = pd.DataFrame.from_dict(model_hash, orient='index')\n",
    "# df.to_csv('../postmodel_derivatives/model_hash.csv', header=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rois = []\n",
    "for model in models:\n",
    "    df = pd.read_csv(datapath + model) #get rois in all tables (some may have failed)\n",
    "    rois += list(df.roi.unique())\n",
    "\n",
    "    \n",
    "rois = list(set(rois))\n",
    "\n",
    "#get inferred\n",
    "theta = df.columns[2:] \n",
    "ntheta = len(theta)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get rois\n",
    "\n",
    "roi_us = np.sort([i for i in rois if i[:2]=='US'])[::-1]\n",
    "roi_other = np.sort([i for i in rois if i[:2]!='US'])[::-1]\n",
    "rois = list(roi_us) + list(roi_other)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.350131004943883\n",
      "0.06874755354613585\n",
      "0.004267655670099832\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJEAAAFTCAYAAACJVp2qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZgUxf3H8c8XdrlvFRABEQRBvAUVBFmCRxC80WhENGrUeERNPDDxF4+gUXN4RqLGqAlGPBAUEI2irKKgIIeC3AhyL4fcN9Tvj+qB3dmZ6ZnZWWbZfb+ep5/e6a6qqdnpp6vn21XV5pwTAAAAAAAAkEilbFcAAAAAAAAAZR9BJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAlpM7OXzczFWDaY2XQze9bM2mW7ngAAAAAAoOQIIiETdkhaESwFkmpIOlLSryRNMbOLs1g3AEAZZ2a/DG5CPJDtumSambWI3GTJdl3SZWa1zGy1mS0ys2rZrg8AlISZPRScl3+R7bpkmpnlBZ9tQbbrki4zO9TMdpjZhGzXBbERREImfOGcaxwsjSRVk9RT0gJJVSS9ZGYHZbOCALC/MrMaZvYrMxtuZj+Y2WYz22Rm35vZW2bW18yqJ1HOE4V6jP4+ifT3x+ltusnM5pjZK2Z2UgY+X3VJ90naIOnJkpaHcGZ2sZl9HASGNpvZDDMbYGa1Y6V3zm2U9JSkppJu3qeVBVDmFBqNMCbGvnhtR/TyRJwyS3WEg5kdLOk2+d8p/ylpeUjMzCqZ2XVmNs7M1gbf52Qzu9PMqsTK45xbKGmQpA5m1mff1hjJIIiEjHPO7XDOvS/p8mBTTUkXZbFKALBfMrNzJM2T9Kyk3pKaSdotaZekFvLn1v9ImmtmP0lQTq6knxfa1C+FauzW3t6mK+RvDhwelDHOzG5LoaxYfi3pEEnPOufWlLAshDCz5yW9Iam7pDryx1JbSb+X7z3cJE7Wp+QDffeYWZ19UVcA+7XotiN6WR8nX2mPcPhDUOajzrmdJSwLCQTXHsMlPSfpFEnVJVWWdJykxySNNbNacbI/LH8MDTCzyvugukgBQSSUpnGSNgZ/H5nNigDA/sbMrpI0TFJjSbMkXSHpQOdcLedcHUn1JPWRNEZSE0mnJSiup6SDJOVLmi2pjZmdkmRVFhXqbdpYvrfpqZKmyF9H/NXMjkrx40mSggvDXwcvX0inDCTPzH4l6ZfyF+Z3SqrlnKst/30ulNRSPsBUjHPuR0lDJDWQdOU+qTCA/VmRtiPG8oc4+UpthIOZ1Zd0laQt8j1dULoGSDpb0lb5/3sN+c4F50haI6mjfICpGOfcHEmfSjoiKANlCEEklDYL1kSQASBJZnaspH/It9PvSTreOTfIObc6ksY5t845N8Q5113SpfK9ROKJ/Oj/r6RXo7alxDm3yzn3haTz5e8YV5LUN52y5C8Mm0ia4Jybl2YZSIKZVZV0f/DySefcX5xz2yQp+D4vkOQknRr0gIvltWB9TWnWFQAiMjzCoa98UGpEMEwXpcTMGku6NXh5t3PuleD6wTnnRki6Oth3mZkdE6cY2pwyiiASSlNn+RO9JM3PZkUAYD8zQFJVSUsk/dw5tyVRYufc65L+FmufmTWQHwq3XdKb8oEkSfpZEFhISzBnwezgZbq9TSOTmhbr/WJmXYP5MApi7KsUzK3gzGxGjP21gkk5nZm1iLH/IDP7k5l9a2Ybg3mepgWTrTZIVGEzO8rM/hXMSbU1qMfnZnZD0HU/JWZ2rJmtCOo6yMxyUi0jSadLaigfKPpr9E7n3GRJHwUvL4/eHxgtabWkY83s+NKoJADEkYkRDonanCuC83CxyZzN7EAz2x3sHxVj/xHBvq0W4+ED5h+w8LSZzQrmodtgZl+b2d1mVjM6fVTeLmY22MwWm9m2YC67j8zsMjOzRHnjlHd60O45M3sk1fwpuEj+OmadpOejdzrn3pG/hjAVHW5f2BD5nrO9zKxhKdUTaSCIhIwzs1wzO0t7u4nukPR6FqsEAPsNMztEUq/g5VPOuXXJ5HPOxXv612XyQwDed8796JybK+lLSfUlnVvS6gbrlHubmlklST2Cl5/HSPKVfBf4g2JMpnqcpLrB321jXFx2lpQj6Qfn3IKo9+0iaaak/pKOkpQrfz3UXtLv5OfcOCJOnW+WNFX+h0gL+fatVvB+AyX9z8xqxP3QxcvrLD8csWGQ/4pSnKOje7Ce5pxbEifNB8E65vxazrld8seOJJ2ZwboBQDJK0uYcKCkS/I7V5uQH6+Ot+EMGTiv03p1jzNHTLVh/6ZzbGvW+F0qaIf9Qgjbygfyqkk6Q9Ij83IKN4tT5UUmfSfqZ/NyBW+Xb7h7yN4T+G7SlSTGzCySNkL/Jf49zrn+yedMQaXM+jf6fFPK/YB2vzVkt317nFCoPZQBBJGRCZzNbHiwr5E9w78tfYO+WdL1zbnE2KwgA+5E87b1YfTcD5UWGrb1aaFuJhrRJ/s6qpNbBy3R6mx6tvRM7T4neGQy1igQsukXtjryODOGLng8qsj+/8EYzO1R+ks8G8kGb1vITfdYM6vM/+cnL347+kWBm50t6WtImSXdJOiiYT6iGpJ9KmiP/3T0e/yMXKe9MSR/Kz231qHPuxgSBwEyI3LmfniDNd8H6oOAHVywTg3XXjNQKAJJT0hEOpwbrZc65ZdE7nXM/yM8NV7lQ2ojCbU4d7Q1GRe+PbnM6ShosHwR5SFJT51xN+Xans/z59GhJ/46uj5ndKt/WrJB0naR6zrm68v+DSyUtD9Z3J/rQhcrrJ98buYqkG51zpdkLSUqtzWmXoFcVbU4ZRBAJmZArqVGwNNTe42qNpJOdcy/FymRmVc3sETNbYmZbzOyroAcTAFRkkV432+Qn1E5b0IOno/yF7/BCu16XD96cFe8OaIIyK5tZJ0lD5c//UnoTlJ4UrOcmGK4XuSCPF0R6OmR/ftT2h+SDNo8EQZu5zrndwTJNfrLPb+Qvfi+IZAoCSpHHUV/snPuzc26VJDnntjvnPpCf+HWzpKvNP0I6ruDO9HD5AFTcu8GW3GOyYy4xiovUaWmCqhXeF+8zTA3WJycoBwCaFbrJHL18FJ7dy+AIh0ib802CNGFtzjMh+6PbnMfl28mbnXP3RnqBBnMDjZN0lqRlks40sw6RTGZWT35Y+1ZJZznnXoj0SnbObQmGsF8o36vpTjOrkuAzycxukfRykL6fc25gjDQtStDmLIjxtqm0ObWCJRbanDKIIBIyId85Z845k5+s7jhJb8nf6X3R/JMQYnlZ0m/lJ027Vb5RGGlm0SdmAKhIDgjWP2agZ0qkp9HQwoEa51yBfC+YHMWf/yaiyA8B+afafCF/rpek+51zX8bPHlfkAnNVgjSfBus97UJwt7KrfGDsSfmL4sL7q8sHzqRCF/TBMLOL5XvIxpw/yjm3Xb79kqQzCu3Kk3So/FCwD6LzBXnnSRov/z/Ni/eBzOwX8vNx5Cr8bnCix2OHLdEid/ATza+1udDf8S7oI9/XgVZ68zcB2P9V0t6bzNFLvJ6OUumNcEi3zakv31tohnxPnuj9reSHmu2Qn7ep8PZTJa2V9GKsN3POrZEUmWOpcJtzkfw5+CPn3NRiGX3ecZK+lx/edmK8D2Rm/yfpKfl5Efs45+Ld9Nml9NublTHKy3Sbk/DmDPYtGn9kVDD8YKqZXSJ/UjxL/tGNlxROZ2YnyXfB7O+cezTY9m9J0yT9WXvvFgAA0hD0nrkiePnfGElelR+GdaXiBFUCkR8C0bZKusg5916aVYz8iPgxQZpx8hfmB5tZ6+CRv0fL36R43zlXYGbTJB1lZgcE8yd0ku+uvzSY/ynixGC7k/RtgvlIqwfrZoW2dQ7WrYNAWjyReZqaxdppZrfJ/693yd8NTtiDyznXONH+LCn8fR0oP6QCAKItdM61SCNfZIRDtDXyvXImxtiXjGTanMiNhw5mVsM5t1n+pkWlYN9U+Ymiu5pZJefcbu0NKE0I0kdE2o1akhYnaHMiwZNYbc5PQtqcyIMgmqlQACtgZvY3SbfLD8M+zzk3Ol5BzrlFkspym5Mo8Ih9jJ5IKBXB3fNfy18oXxyjd1Ef+bsJzxfKs1U+Ut/RYjxNBwAqiNXBun6COQKScbqkJpIKtPepW4UNk78LeIyZHZugnIWFeptWkdRWfj6hapKeK8H5OvJkuO3xEgQX5JEn5XSLWo8J1vnyc0h1jdofuaMcEbmLaYp/d7yR/HwXkh9qFp23akjeajHyFvZ48P4PhgWQSsGmYF09QZrC9Y73+OvCE6QmKgsA0pHuCIcwybQ5c+WHWOVqbyBnT5sTBI0+kx8WfUzU/nhtTo4StxuRHjux2pwaIXlzY+SNaC4fQJKkXyUKIJWSTLc5VVKZRByliy8CpcY5N1t7xyw/FLX7eEnznHPRdwO+KrQfACqiyCPrq0qK+ZSwJEWGsjWUtDPGnDkbtPcCLqkJtp1zO5xzs5xzN0p6QVJTSa+leWG3JljXC0kXPbwgeu6J6Dks4s1NEanjusgPlJAlL0bed5LMe3+czzI4WN8R9MjdlyJzTzRJkKbwvmITzwYK/4BbHScNAJSYc25bMJzrEvmnRx4jP8IhHdlqc6Ym2W5cFSPvk0nmfTnG51heqE5/CobX7UuptDkbnXMb4qSJtDlrgyAeygCCSChtfwnWp5pZXqHtByv2BWpkW6ITDgCUZ/nyQ64k6dx0CjCzOpLOTyHL5WnMb3O3fLf+U7R32FwqIvMchN1Vjr5gP03+jmVkSMOeC34zq6q9k29GX9BH5gmqY2Z1lZpI3uYp5ot2haS35Xs7fWBmCW+YJJiUNnSJUVzkKTjtE7xl5Gk6KyMTh8cQ+b62O+fWJ6o/AGRCEiMckpFymxO0pcdJmuWcWx5j/6Hy8+XtkvR5VDmRdiPm8OYQmWhztknqLV+vQyR9HNQ3JjNLNBF62DIhRpGptDkzEqSJfF+J5rLCPkYQCaXKOTdZe4dR3FtoV3X5k1u0rYX2A0CFE0wYGpln6JbgIjZU1NC3S+TPo4vkL8DiLQfI9yZpKD8/Uir1/FHS34OX96cRhIo8ea5FSLrP5S/Qm5lZb0kHSfrCObczqEeBpJmSjpWfh6+apALnXPRF6URJO+WHk6X0WbV3roljzOyQFPPuEdT5Uvkns9WT9KGZHZ0gS6JhDGFLtE+CdXuL//S4M4N1omEPLYJ1iZ4cCACpCBnhkIzIOeuwkHSRGxMnyU92XVlFb0pMkr+RcZr2PkRhUoyeNJF2o4GZpfpksUjevOBhEWlxzm2UdLb8SI/m8oGkpnGSV1b67c1BMcqLtDldzaxajP3S3snEk2lzZiZIg32MIBL2hceCdQ8zOyX4e4v2jk0urFqh/QBQUd0rH2hvKum/CS7AJEnBwwx+U2hTZHja2865tQmWNZLeicqTiqeDeraQ1DfFvF/I97iqn6ibfXBhPjl4+YdgPSYqWb78NU3kZkX03BSRcoYELx80s9rx3tPMcsys8JNiRssH5CrLP/whrrD5OpxzO+SfEve+fBDvIzNrFydtMsMYYi4xihstPz9WJfkno0bX+1j5ebQkP+l6PJEn341N9DkBoBTEG+GQjEhPofaJ2lTn3HfyTxurKt/jVirU5jjnIr2ODpB0U7A5uuernHMz5Z/YKUmPmVludJoIM6se9KSNeFN+TqH62tvuxcsb1uasl7/BMklSS/lAUrEbCc65BSVoc1rEeOu35a8P6km6Nka9z5Efsu/kn9QdD21OGUQQCaXOOfeh9v4A+L9gvUyxH9UY2bY0xj4AqBCcc1PkL06dpF6SJptZXzOLPIlFZlbXzC40s0/k787WDra3ktQlSPZ2Em8XSXNOqhOWBt37/xO8vCeVuZGCAFakt1DHRGm1NygUSRd9wZ4fsj+iv/y8GG0kfWFmP41c2JvX2sx+I3/Hs0Ohuu6QdLP893GZmQ0zs+Mi+80s18w6mNlj8o9cTsj5J5leIB/YaShptJm1DstXEsF73h+8vN3Mfhv50WJmnSQNlb8u/Nw5NyJBUZH/cbFAHQCUpgQjHJLxjaT18pNRHxeS9rNgXdI259fygZTT5M/zXSLtpJlVNrOjzewPkuar0O8i5580ek/wsr+ZvWBmbSL7g6BTVzMbKH9DJiHn3Fr5Xj/fSGod1CVW76GMCa4PngxePmZmV5h/aqzM7GxJLwX7XnPOfROrjKCH9YnBS9qcMoQgEtLmnLsqxuSj8dKeEKTtFWyaIqlVjB8sJxfaDwAVlnPuRUkXyvceaSsfrFltZhvMbL2ktfI9a/IkLZT0cZC1X7BeoeTu3H0of2FdVX6oVar+Iv+0zTaSfpZi3sjQhF4JUxW9QC/8xLZY+2O9luTvtMoPZVsq6ShJoyRtMrNV8sOpZ0v6q6RW2jsvVSTvu5KukX+yz3nygb3NZrZavvfsBEl3SkpqviXnn0h6blDXg+XvDocNsygR59xA+QnRK8l/bxvMbIP8j5DD5H/IXBIvv5k1k3S0/B3ykaVZVwCII9YIh1BBD6JIb9RU2pw5zrnom9uF9+9WnLbWOTdB/obBOvkniH4maXPQ5myRD+o8IKmxirc5T8vffHfyPXlmmdlGM1sjP5zuU0k3aO8ojoSCGzeny89V1E6+F+wByeQtgXvlh+dXl/Rv+fY20n4cIN9u3pAg/6nyPZkWSPqyVGuKlBBEQra8JX/8XRfZENwR/YWkr51zoXdyAaC8c84Nk+9+fpP8hdhi+ccF58hfVL0l6eeSjnDOfRrctYsEkd5J5kkmzrntkiI9T1Ie0uacmyXp3eDl76LmZgrzkvwF+HkhQ/Y+C9JJfj6kHVF1WCppbvByjaRpCeo7QT4od7d88GSj/EXqZvl5k56S1M05F2t4wkvy3e+fkDRdfq6mOvLzSo2RdJ9SeKKec26z9k582lTSJ2ZW0sm7w97zOvlg3yfynz1HvufVQ5KOi/FjqbCfyc8p9UaCJ+kAQKmJM8IhWS8G67AbHvlx/o6YIN9mSNI3QU+fmJxzo+RvsgyQH1IWGeK1Xr4NekTSic65hTHyDpCf7+95SXPkfzvVlB/R8YGku+SDU0lxzq2U1EN+fqhj5OflC3taXdqCtvoc+UDRePnP7uQ7C9wtqUtIWxK5sfVSMLk6ygjj+0C2mNkb8tH5J+RPjP3kn/JzhnPuk0R5AQDlg5mNkL8rfLFz7q1s1wfxmdnXkk6Q1Nk5Ny4sPQCUNWY2Tf6JYR2dcxPD0iM7god1LJbvsdTSObcoy1VCIfREQjb1k/S4pMvl7/xWk3QOASQAqFDul78zWWyyZ5QdZtZdPoD0AQEkAPux+4L1HVmtBcJcIf/ktxcJIJU99EQCAABZZWavy8/Fc4Zz7qOw9Nj3zGy0pO7ywy4mh6UHgLLKzMbLT4rdzjk3O9v1QVHB5OMzJDWTdHjIMGtkQU62KwAAACq8/vIXjLWyXREUZ2a15CdxHUwACUA58Cv5ByQ0kX+gAsqWJpJekzSdAFLZRE8kAAAAAAAAhCpJTySiT+Xcz57zUx68fn2nLNcEFUwqT3YCaIuyhDYC5RxtEZJFO1QO0cahDCiz7RATawMAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQKifbFUBqBg0ekfEy+17aO+NlAgCKGjRmZLarkFEr1vp19Ofqm9crC7UBAJRVYycMyXYVUrZuQxVJ8evepeNF+7I6QJlCTyQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAiVk+0KIPsGDR4Rc/uKgsT74+l7ae+SVgkAAAAAAJQx9EQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACBUTrYrgPJn0OARGS+z76W9M14mAAAAAABIHj2RAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACBUTrYrgP3XG4P/reHvvClJuvTyX6hX7wuTzvu/94dr1szpWrxoodavX6ctWzarRo2aan7oYep6Wg917pInMyuWLy8vT/n5+XHLPeuss/T++++n/mEAAAAAAEBCBJGQlvnzZmvk8CEyMznnUs4/YvgQrV+3Tk2bNdfhbdqqatVqWr2qQN9N/0bTp03VV199rltv/50qVYrdWe6ss85S48aNi20/+uijU64LAAAAAAAIRxAJKduxY4eeG/iE6tatp5at2ujrieNTLuOmW+7UoS1aqVq1akW2L160UI88dK8mTfxSYz/9WKflnR4zf//+/ZWXl5dO9QEAAAAAQBqYEwkpG/Lmq1q6ZJGuuuYmVa9RI60yjmjbvlgASZKaNjtUp5/ZS5I07dvJJaonAAAAAADIHIJISMncubM0auRQdTq1m0448aRSeY9KlSpLknJyc0ulfAAAAAAAkDqGsyFp27dv1/MDH1fNWrV1Rb9flsp7FBQs18ejR0mSTjjx5Ljphg4dqqFDh2rbtm1q0qSJunfvrq5du5ZKnQAAAAAAAEEkpOCt1/+jZUuX6KZf36XadepmpMxPx3ykGTO+1a5du7Rm9SrNnTNTu3c7nXv+xerQsVPcfE899VSR1/fdd59OPfVUvfbaa2rWrFlG6gYASN2gMSOzXYVS1zevV7arAAAAkBUEkZCU2bNn6P1R7+rEDqfolE6Z6/Eze/Z3Gvvpx3teV65cWX0uuVw/Pfv8mOm7du2qfv36qWvXrmratKlWrlypL774Qr/73e/0+eef6/TTT9ekSZNUs2bNjNURAAAA2BfGThiS7SoAQEIEkRBq+/ZtemHgE6peo7quuvpXGS372ut+rWuv+7W2b9+mlQUr9Gn+R3r7rdf05fixuuOu+1S/wQFF0v/xj38s8rp58+Zq3ry5evbsqRNOOEGzZ8/WwIEDdccdd2S0ngAAAAAAVHRMrI1Qbwz+t5YvX6rL+16revUblMp7VKlSVYc0ba7LLr9al1zaTz8s/F7/fvm5pPPXrVtXt956qyTpvffeK5U6AgAAAABQkdETCaG+njBeZpX02aej9dmno4vsW7Z0sSTp4w9HacqkCWrU+GBde92vS/R+Xbv10Guv/kuTJ32lnTt3KicnucO0bdu2kqQlS5aU6P0BAAAAAEBxBJGQFOd2a+aMaXH3FxQsV0HBcm3evKnE71WzZi1VrlxZu3bt0qaNG1S3Xv2k8q1evVqSVKtWrRLXAQAAAAAAFEUQCaEef/rFuPueG/i4xn76sS69/Bfq1fvCjLzfrJnTtWvXLtWoWVO169RJOt8bb7whSerYsWNG6gEAAAAAAPZiTiSUmtdfe0V3/fYGvf7aK0W2z5o5XZMnfaVdu3YVyzN71nf653NPSZK65Z2hSpUq79k3ZswY5efnyzlXJM/mzZt11113adiwYcrJydEtt9xSCp8GAAAAAICKjZ5IKDVr167RsqVLtHbtmiLbV6xYphf+8aRq1KypFi1aqW69+tq6ZYsKVizTkiWLJEnHHd9BfS7pWyTflClTdPvtt+vggw/WscceqwYNGmjFihWaMmWKVq9erapVq+rFF19U+/bt99lnBAAAAACgoiCIhH2ubbujdP6FP9Osmd9pxfKlmjN7piSnunXrq+NJndW5S546dOxULF+3bt10ww03aOLEiZo8ebLWrFmj3NxctWjRQpdddpluueUWtWnTZt9/IAAAAAAAKgCLHhqUgrQzIn2DBo/YZ+/1zzl+fW3rffaW+0zfS3tnuwqIz7JdAexX9pu2aNCYkdmuQkb980u/vvbk7NYjG/rm9cp2FVD6aIuQrIy2Q2MnDMlkcUjTgPwqkqR7u22Pub9Lx4v2ZXVQMZXZdog5kQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRJgiN4AAB8QSURBVBAJAAAAAAAAoQgiAQAAAAAAIFROtitQ3g0aPCLbVQAAAAAAACgxeiIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABC5WS7AmXJoMEjsl0FAAAAAACAMomeSAAAAAAAAAhFEAkAAAAAgBAfvv+Zbvrlvapbt65q1aqlDh066O9//7t2796dVnnvv/++zjzzTDVo0EA1atTQUUcdpYceekjbtm1LmO/LL7/UBRdcoIYNG6patWpq3bq17rrrLq1bty7p9542bZqqVq0qM9NRRx2VVv1RMRFEAgAAAAAggb899oIe/MOTmjlzvrp27aozzjhDs2fP1s0336w+ffqkHEh67LHH1LNnT3388cc64YQT1KtXLxUUFOjee+9VXl6eNm/eHDPfa6+9plNPPVXDhg1TmzZtdN5552n79u3685//rA4dOqigoCD0vXfu3Kkrr7xSO3bsSKnOgEQQCQAAAACAuMZ8PF5D3/pADQ6op5df/atGjBihoUOHas6cOWrXrp2GDh2qp59+OunyJk6cqP79+6tGjRr6/PPP9dFHH+nNN9/U/Pnzddppp2n8+PH6/e9/Xyzf4sWLdc0118g5p2HDhmns2LF6/fXXNW/ePP3sZz/T3Llzdf3114e+/8MPP6xJkybpxhtvTOn/AEgEkQAAAAAAiGvQy29Lkn51c181a37wnu2NGjXSwIEDJUmPPPJI0r2RHnnkETnndPfdd+vkk0/es71WrVp66aWXVKlSJT377LNau3ZtkXxPPPGEtmzZoiuvvFLnnXfenu05OTl6/vnnVadOHQ0bNkzfffdd3PeeOnWqBgwYoAsvvFB9+vRJqr5AYQSRAAAAAACIoWDFas2aOV+5uTnq3qNTsf3dunXTIYccouXLl2v8+PGh5W3fvl2jRo2SJF1++eXF9rds2VKdOnXS9u3b9d577xXZN2zYsLj56tSpo3POOadIumg7duzQVVddpdq1a+vZZ58NrSsQC0EkAAAAAABimDP7e0lSi5bNVLVa1ZhpOnbsKEmaPHlyaHmzZs3S5s2b1aBBA7Vq1Srp8tavX6958+YV2Z9qPQYMGKApU6bo8ccfV6NGjULrCsRCEAkAAAAAgBiWLV0hSWrc+KC4aZo3by5J+v7770PLi6SJ5Em2vAULFkiS6tWrpzp16qRcj8mTJ+vhhx9Wz5491a9fv9B6AvEQRAIAAAAAIIbNm7dKkqpVj90LSfJzGUnShg0bQsvbuHGjJKlmzZoplZduPskPobvyyitVvXp1Pffcc6F1BBLJyXYFAAAAAABA6XjwwQf17bffauDAgWrWrFm2q4P9HD2RAAAAAACIoUaNapKkrVu2xU0T6SVUu3bt0PIivYU2bdqUUnnp5vv666/16KOPKi8vT9dff31o/YAw9EQCAAAAACCGxgc3lCQtX74ybppFixZJklq0aBFaXiTNDz/8kFJ5hx56qCRp7dq1Wr9+fcx5kWLlGz58uHbu3KkVK1aoe/fuRdKvXbtWkp9DKS8vT5L0z3/+U4cffnjo50DFRRAJAAAAAIAYWh9xmCRpwfxF2rZ1W8wntE2YMEGSdPzxx4eW17ZtW1WvXl1r1qzRvHnzYj6h7auvvipWXt26ddWqVSvNmzdPEyZMUI8ePZLKFzFjxgzNmDEjZp02b96s/Px8SXt7MwHxMJwNAAAAAIAYGjU6UG3attSOHTv1yehxxfbn5+dr8eLFaty4sTp16hRaXpUqVdSzZ09J0quvvlps//z58zVu3DhVqVJFvXr1KrLvvPPOi5tv/fr1Gj58uCTpggsu2LP9/vvvl3Mu5vLJJ59Iktq3b79n23HHHRf6GVCxEUQCAAAAACCOvlf6oMzAZwZp8aJle7YXFBToxhtvlCT1799flSrt/Xn9zDPPqG3bturXr1+x8vr37y8z06OPPrqn95DkewFdffXV2r17t2688UbVq1evSL7bbrtN1atX1yuvvKJ33313z/adO3fq+uuv1/r163X++efryCOPzMwHB2JgOBsAAAAAAHF079FJ5190loYN+UBX/vy3evWMD5Sbm6vRo0fvCdzcfPPNRfKsWrVKs2bNUuPGjYuV17FjRz3yyCO6++671blzZ/3kJz9RvXr1lJ+fr4KCAp188sl66KGHiuVr1qyZXnzxRV1xxRU6//zz1aVLFzVp0kTjx4/XwoULdfjhh+u5554rtf8DINETCQAAAACAhH579y/1hwdvVZsjDlN+fr4++OADHX744XrmmWc0ZMgQVa5cOaXy7rrrLo0aNUrdu3fXhAkTNHz4cB144IEaMGCA8vPzVaNGjZj5LrvsMn3++ec699xzNWPGDA0dOlQ5OTm68847NXHiRDVs2DATHxeIy5xz6eZNO2NZNWjwiGxXoUz55xy/vrZ1dutRGvpe2jvbVUB8lu0KYL+y37RFg8aMzHYVMuqfX/r1tSdntx7Z0DevV3gi7O9oi5CsjLZDYycMyWRxSNOA/CqSpHu7bY+5v0vHi/ZldVAxldl2iJ5IAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUT2cDAAAAACBJ5XHuKuZ5QrLoiQQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhMrJdgWAbBg0eERGy+t7ae+MlgcAAAAAQFlDTyQAAAAAAACEIogEAAAAAACAUAxnAwAAGfHFR2M0+t33tGj+Au3etVtNmjdV15+erh7nna1KlVK/b/XNV19r1JvD9P2sOdqxfYcOOriROvXoprMvuVC5VXIT5h33cb4+G/WRFsydry2bNqlWnTo6pEVzdT49T6f99PRi6Teu36CRrw/RpLHjtXJ5gXKr5KpZyxbK63WWupz5k5TrDgAAUB4RRAIAACX28hMDNfqdkcqtUkXtTzhWlXMqa/qkqfr3U//Qd5On6pb770kpkDTitbf0+vMvq1KlSmp33NGqUbuWZk6dprde/I+mjPtK/f/6kKpWq1Ys3/bt2/X0fX/SlPETVLVaNbU5qp1q1qmtH1et1vyZsyXnigWRCpYu159+8zutWlGguvXr6egOx2vzpk2aN2O2Zn0zXdMnTdV1d98mMyvx/wkAgLJo7IQh2a5CxnXpeFG2q1AuEUQCAAAlMiH/c41+Z6TqNqive598RI2bHiJJWrfmRz38m99p4mfj9OHbw3VWn/OSKm/+rDl644VXVKVaVd3z14d1+JFHSJK2btmiv/R/QLO+maY3X/yP+t70y2J5n3/kcU0ZP0End++qX/zmJtWsVWvPvh3bd2jJgoXF8vz9j49p1YoCdex2qq7vf/ue4NSShYv0l7vv09gPRqvNUe3UvfdPU/7fAAAAlCfMiQQAAEpk+H/flCRdet1VewJIklS3QX1ddduNPs1rb2n37t1JlTfiv2/KOafel/bZE0CSpGrVq/seQZUqafQ7I7Vp48Yi+b756mt9+clnat7qMN34+zuKBJAkKbdKrlq0ObzItjnTZ2j+zNmqUbOmrvntLUV6Nx1yaDNd9qurJUnv/Od1OeeSqj8AAEB5RRAJAACkbc3KVfp+9lzl5ObopLwuxfa3O+5o1T/wAK1b86PmfjcrtLydO3Zo6ldfS5I6n55XbH/DJo3V+si22rljp6aOn1hk34fDRkiSzrroXFWqXDmp+s+fOUeS1KJNK9WsXavY/qM7nCBJWl2wUvNmzE6qTAAAgPKKIBIAAEjbgjnzJEmHtDhUVapWjZmmZdvWkqSFQdpEli1aou1bt6lWndpqdMjBMdMcFilv7vw923bv2qXvJn8jSTrimPZas3KVRg4eopf+9oz+O/BFTcj/XLt27SpW1rYtWyVJtevWjfle1WpUV06uH/2/YPbc0PoDAACUZ8yJBAAA0rZy2QpJ0oGNDoqb5oCGft/K5SuSLi+SJ2F5y5bv2bZi6XJt37pNkjTr2+/0ypMD97yWpFGSmhzaTL956P/U6JAme7bXqe+DRwWFyipszcpV2rljZ9L1BwAAKM/oiQQAANIW6ckT60lpEdWqV5ckbd28JbS8rVu2JFFetWLlbdqwYc/f//rrM2rdvp3++PyTeuG9N/XAwL+pdft2Wrpwkf5yzwPasX3HnrTtjjtGZqYFs+dq/qw5xd5r9Dvv7fl7y6bNofUHAAAozwgiAQCA/d7u3XsnvT6g4YG640/3qUXrVqpWvbpatm2ju/78oOo2qK/li5Zo3Ogxe9I2OuRgdT49T845PXHvHzXxsy+0acNGrVpRoGH/fk3vvf62Kuf4jttWyfb1xwIAAChTGM4GAADSVjXoFbRt69a4aSK9i6rVqB5aXqTXUuLythYrr3qhv7uc2UM5ubnFyj31jO567/W39d3kb3RazzP27Lvq9pu0dcsWfT12vJ78w8NF8p2c11U7d+7Q12PHq1bt2qH1BwAAKM8IIgEAgLQd1LiRJGnVipVx06xZuSpI2zC0vAODNKsLEpQX7DsweO/C+STpoIMbFctTuK7rflxbZHu16tV02x/v1ZzpM/TNV5O0dvUa1apTW0d3PEFHHn+MHrj5DklS05aHhtYfAACgPCOIBAAA0nZo65aSpCULFmr7tm0xn9A2f6afa+jQw1uFltekeVNVqVpVG9dv0Ioly2I+oW1veS33bKteo4YaN22i5YuXauP6DcXySNKGdesl7Z1TKVrr9u3Uun27Itu2bN6sH+Z+r8qVK+vI444JrT8AAEB5xpxIAAAgbQc0PEgtWrfSzh079dWYscX2z5jyrdasXKW6Derr8PZtQ8vLyc3VMSedKEn64qMxxfYXLF2uOd/NVE5ujo7r1LHIvg5dO0uSpk+aErPs6ZOmSpIOa9M6tB4Ro995T9u3bdNJ3bqoboP6SecDAAAojwgiAQCAEjnn8oslSYOff1krlizds33dj2v1yhMDfZrL+qhSpb2XHR8OHa67+t2gfzz81+Ll/byPzEwjBr+leTNm7dm+dcsWvfDYE3K7d6vHeb1Us1atIvnOuuhcVateXVPGTdCnoz4ssm/Um8M065tpqlqtmk7reXqRfct+WKxNGzYW2eac05iR/9Nb/xqkWnVq6+c3XpPKvwQAAKBcYjgbAAAokZO6dVGPc8/W6Hff0z1X36z2Jx6rnJwcTZ80VVs2bdaJXU7RGRf0LpJnw7r1WrZoseo2qFesvJZt2+iSX16p159/WQ/efKeOPOFY1ahZUzO/mab1P65Vq3ZH6OJrriiWr94BDXT9PbfrmQce1QuPPakPhryrRk2baOmCRVqy8Afl5ubqht/9RvUOaFAk37iP8/Xuq2/qsDat1OCgg7R79259P3uOVq9Yqbr16+mORx8olgcAAKAiIogEAABK7Krbb1Sbo4/UR8NGaubUadq9e7eaNGuq03qeoR7nnV2kF1Iyel/WR81bHab33hiq+TNna8f2HWp4cGOdeeE5OvuSC5VbJTdmvg5dO+uBfzyu4a++qZlTv9WShYtUu24dderRTef8/GI1a9miWJ4jjz9Wi79fqO9nz9MP8xaoUqVKOujgRjrtytP10z7nq0atmun8SwAAAModgkgAACAjOp+ep86n5yWV9sKrLteFV12eMM0xJ524Z36kVBx6eEvdfN/dSadve+xRanvsUSm/DwAAQEXDnEgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAINR+PSfSoMEjsl0FAAAAAACACoGeSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAg1H49JxIAAAAAAEC0sROGZLsKaevS8aJsVyEueiIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFAEkQAAAAAAABCKIBIAAAAAAABCEUQCAAAAAABAKIJIAAAAAAAACEUQCQAAAAAAAKEIIgEAAAAAACAUQSQAAAAAAACEIogEAAAAAACAUASRAAAAAAAAEIogEgAAAAAAAEIRRAIAAAAAAEAogkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIFROtisAILZBg0dktLy+l/bOaHkAAAAAgIqFnkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIRiYm0AAIAUDBozMttV2Cf65vXKdhUAAEAZQ08kAAAAAAAAhCKIBAAAAAAAgFAMZwMyYNDgEdmuQlZk+nP3vbR3RssDAAAAAGQOPZEAAAAAAAAQiiASAAAAAAAAQjGcDQAAAMXwFDoAABCNnkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQBJEAAAAAAAAQiiASAAAAAAAAQhFEAgAAAAAAQCiCSAAAAAAAAAhFEAkAAAAAAAChCCIBAAAAAAAgFEEkAAAAAAAAhCKIBAAAAAAAgFA52a4AAAAAkC2Dxowstq1vXq8s1AQAgLKPnkgAAAAAAAAIRRAJAAAAAAAAoQgiAQAAAAAAIBRBJAAAAAAAAIQiiAQAAAAAAIBQ5pxLL6PZdc655zNcH5RTHC9IFscKUsHxgmRxrCAVHC9IFscKUsHxgmSV5WOlJD2RrstYLVARcLwgWRwrSAXHC5LFsYJUcLwgWRwrSAXHC5JVZo8VhrMBAAAAAAAgFEEkAAAAAAAAhCpJEKlMjs9DmcXxgmRxrCAVHC9IFscKUsHxgmRxrCAVHC9IVpk9VtKeWBsAAAAAAAAVB8PZAAAAAAAAECrlIJKZ/dzMPjOzdWa20cwmmtlNZkZACnuY2ctm5hIsM7NdR+w7ZnaEmd1qZoPMbKaZ7Q6Ogz5J5OWcU0Fk+rs2s5+a2f/MbI2ZbTazaWb2ezOrmum6Y9/K1LFiZs3M7Fdm9qKZfWNmO4Nz0x2lVXfse5k4Xsyskpl1NrMBZvaFmf1oZjvMbIWZvWdm55fmZ0Bqykp7YmYnm9lQMysws61mNsfMHjOzuiH5jgiumZaa2TYzW2hmA83s4HTqj8T21+PFzPJCfm85Mzslnc+A2LJ9rJjZgWZ2dXA+mBCcH5yZPZPk+6V1TipWTirD2czs75JulLRV0mhJOyT1kFRb0lBJfZxzu1OpAMonM3tZ0pWSPpc0N0aSZc65e/ZppZA1ZvaEpFtj7LrYOfdWgnyccyqITH/XZnaXpEcl7ZI0RtKPkrpJOkjSeEk9nHObM/gRsI9k8lgxs9skPR5j153Oub9kpsbIpkwdL2Z2uKQ5wcs1kibKn1daSuoYbH9Z0tWOuSKyqqy0J2Z2maT/SKosfz28RNIpkprLXxuf6pwriJGvm6RRkqpLmiR/3B0rqa2klZK6OOdmJ1t/JLY/Hy9mlifpE0krJL0fp0p/dM7NS7b+iK8sHCvBDYuhMYr7u3Pu5pD3S+ucFJNzLqlF0kWSnKRlkloX2t5I0nfBvluTLY+lfC/yF1JO0lXZrgtL9hdJ10p6TNIlkloFJ0onf7KNl4dzTgVZMv1dS+ogabekTZJOLrS9lqT8oLzHs/25WcrEsXKepCckXSGpnaR/B2Xcke3PylK2jpeg7Rot6aeSKkft6yZpY1DeL7L9uSvyUlbaE0lNJW2W/3F4XqHtOZIGB/mGxshXM6i7k3Rz1L6/BNu/VtARgKXCHy95wb4x2f5flvelDB0rnSQ9K+kaScdJGhCkfSbk/dI6xuKWl8IHnRgU3i/Gvm6F/qmVsv0ls2R/EUEklgSLkgsicc6pIEumv2tJbwV5/hBjX8ugAd0mqV62PztLdo+VGGVE2i6CSOVg2ZftiKR7g/JGZ/tzV+SlrLQn2hvw+VeMfHUkrQv2Hxm17+Zg+8cx8lWW7y3gJJ2d7f91eVjKwfGSJ4JIFepYiZH2fiUXRErrGIu3JDV2z8yaSjpR0nZJb0bvd87ly3eHaizfJQoA0sY5p+LI9HdtZlUk9QxevhqjvPmSxkmqIunstCuOfY7zAlKRheNlcrBumoGykIYy1p5E5siKlW+9pOFR6ZLJt0u+x0CsfEhROTlesA+UsWMlXRk9xpKdAOr4YD3dObclTpoJUWkBSepuZn8zs+fN7I9mdla6E4+hQuGcU3Fk+rs+QlINSWtc/DkAOHb2T5wXkIp9fby0DtbLMlAW0lMm2hMzqyM//LHw/mTrcXzU/mTzIXXl4XiJaGRm9wW/tx4PJl4+IIk6Izll4lhJV4aOsSJyknzvw4L1wgRpfohKC0hSvxjbvjOzS51z3+7z2mB/wTmn4sj0dx1J80OCNBw7+yfOC0jFPjtezKyGpF8HL4eUpCyUSFlpT1oE67XBHf6k8gU/9BoEL+N9Bs5xmbNfHy9R2soPayrsaTPr75x7OkF9kJyycqykq0WwLskxVkSyPUJqBetNCdJsDNa1kywT5dsU+QuqI+WPnyaSekuaGmz7yMwOyV71UMZxzqk4Mv1dc+yUX3y3SMW+PF6elb/w/k7S8yUsC+krK+1JSfMlyss5LnP29+NF8vPYPC6pq/xQqtqSTpD0T0nVJD1lZtcmrjaSUFaOlXRl/P2S7YkEpMQ590TUpk2SRprZh/Izzp8i6R75CQQBAAD2O2b2f5KulP8xd4lzbluWqwSggnDOTdbe+dgiJkv6pZl9I+kpSY+a2X84NyGTku2JFIlM1UyQJhLh2pB+dVDeOee2S/pT8JJJbREP55yKI9PfNcdO+cV3i1SU+vFiZr+R9GDwXj2dc9PTKQcZU1bak5LmS5SXc1zm7O/HS5i/S1olP0Ty5BTyobiycqykK+Pvl2wQaUGwPjRBmmZRaYF4ZgZrhrMhngXBmnNO+bcgWGfqu46kaZ6h8lB2LAjWnBeQjAXBulSOFzO7RdJfJW2R1Ns5Ny7VMpBxC4J1ttuTyLwp9YJ5jpLKF8xV8mPwMt5n4ByXOQuC9X55vIRxzu2WNCd4yW+uklkQrLN9rKQr48dYskGkSDe59mZWPU6ajlFpgXgiTwvYmDAVKjLOORVHpr/rmfI/6hqYWas4aU5KoTyUHZwXkIpSO17M7Cb5YSJbJZ0bPN4Z2Vcm2hPn3DpJkScudSyWI06+wKQ08yF15eF4CcNvrswoE8dKukrjGEsqiOScWyR/Uqsi6eLo/WbWTVJTScslcScGYS4J1vEeMYgKjnNOxZHp7zoYMjsqeHl5jPJaSuokabukkWlXHPsc5wWkorSOFzO7QdIzkrZJOt8591FGKowSK2PtyTsJ8tWRdE7wcmgK+SpLujROPqSonBwvcZnZsZLaSHKSJiabD//f3t2DRhGEARh+YwKCBC2EQLoIQlqJKFiZTiuxSKOVjVVAC42NjbUKWohYiQhyEkFFG3+aICppFJuoqIWKYkIQf0CNaHIWM8J63mXvjphZ994HviKzu5Od3S+zd9/ldv9WsFxp19LmWLVabSqAEUISvgPWZ9r7gKm4bH+z/RnlDWAD4Uls3TXtPcABYD7my7bU+2oky5GJmAMji6zjnNMh0c65JtyU/ylwvk5/m4AFwg39N2faezO5dyL1uI30uVKn/3Oxj4Opx2oUL1+AvXFumSPcAyn5GI1/fs7bup4QvhryNb7m3ZFp7wEqcbsrdbbrjfteBUZrlh2L7Q+BrtTHugxRgnzZB6yt076F8FW2KlBJfZzLEEXJlTr9HInrnspZr60ca9hfiwfvdPwF34DrwGXC0yiqhKpVdyv9GeUMYGfMiffAbeACcAN4G9vngbHU+2ksa04MAZOZ+Bxz4Vm2vc52zjkdEq2e68xFc6JBf4fi8p/ALWAcmIltk8Cq1GM20ucK0F8zN83GdV/VtPenHreRNl8IH5AtxGVPCAXHenE89Zg7PYpyPQF2xW0WgDvARcL9RqqEN/h9DbbbSniz9/s/SCrA4/jzLDCY+hiXKf7nfAE+Aj8I3+4YBy4BjzJz1V1gdepjXJYoUK5kX5+8ietP17QPLUWONTwWbRy83cA9wpvAL8ADYBRYkfrEGsUIYB1wErhPKBzNxT+258BZYGPqfTSWPSeG4wS1aDTY1jmnQ6KVc513YY7rbCcUsj/EOWgKOAysTD1Woxi5Agw0MzcBA6nHbKTNl2avY8DL1OM1inM9ITwV6yqh+PMdeAEcBdbkbDdI+BB2Om73GjiDBW3z5c/1x4BrhPvdfCIUlGaAm8Ae/LC1lLnS5LVoeClyrFF0xc4kSZIkSZKkhpp9OpskSZIkSZI6mEUkSZIkSZIk5bKIJEmSJEmSpFwWkSRJkiRJkpTLIpIkSZIkSZJyWUSSJEmSJElSLotIkiRJkiRJymURSZIkSZIkSbksIkmSJEmSJCmXRSRJkiRJkiTl+gVHxNHaA46migAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "font = {'family' : 'sans-serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 22}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def simpleaxis(ax):\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    return\n",
    "             \n",
    "fs=24\n",
    "\n",
    "theta_ = [\"R0\",\"car (week 0)\",\"ifr (week 0)\"]\n",
    "label_ = [r'R$_{0}$','CAR (week=0)','IFR (week=0)']\n",
    "xticks_ = [[0,5,10],[0,0.1,0.2],[0,0.005,0.01]]\n",
    "# xlim_ = [[0,5,10],[0,1,2],[0,0.5,1],[0,2,5]]\n",
    "\n",
    "def histrois(ax,theta,label,histcolor,xticks):\n",
    "    x = []\n",
    "    dfbest = pd.read_csv(\"../postfit_derivatives/\"+theta+\"_summary.csv\")\n",
    "    for roi in rois:\n",
    "        model = dfbest.loc[(dfbest.Region==roi), \"Model\"].values[0] + \"_fit_table.csv\"\n",
    "        df = pd.read_csv(datapath + model)\n",
    "        try:\n",
    "            x2 = df.loc[(df.roi==roi)&(df['quantile']==0.5), theta].values[0]\n",
    "            if np.isfinite(x2):\n",
    "                if theta in [\"car (week 0)\",\"ifr (week 0)\"]:\n",
    "                    x.append(x2)\n",
    "                else:\n",
    "                    x.append(x2)\n",
    "        except:\n",
    "            print()\n",
    "    mu = np.median(x)\n",
    "    print(mu)\n",
    "    f = sns.distplot(x,hist=True,kde=False,ax=ax,color=histcolor)\n",
    "    simpleaxis(ax)\n",
    "    ax.set_title(label,fontsize=fs)\n",
    "    ax.axvline(mu)\n",
    "    ax.text(mu,20,str(np.round(mu,3)))\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "    ax.set_xticks(xticks)\n",
    "    ax.set_xlim((min(xticks),max(xticks)))\n",
    "    return\n",
    "    \n",
    "f,ax = plt.subplots(1,3,figsize=(20,5))\n",
    "c_ = sns.color_palette(\"cubehelix\")\n",
    "for i in range(len(theta_)):\n",
    "    histrois(ax.flatten()[i],theta_[i],label_[i],c_[i],xticks_[i])\n",
    "        \n",
    "plt.subplots_adjust(wspace=0.5)\n",
    "plt.savefig('../postfit_derivatives/fig_regionaverages.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
